/* ************************************************************************* **
**    OpenVFIFE - Open System for Vector Form Instrinsic                     **
**                Finite Element Method (VFIFE)                              **
**                GinkGo(Tan Biao)                                           **
**                                                                           **
**                                                                           **
** (C) Copyright 2021, The GinkGo(Tan Biao). All Rights Reserved.            **
**                                                                           **
** Commercial use of this program without express permission of              **
** GinkGo(Tan Biao), is strictly prohibited.  See                            **
** file 'COPYRIGHT'  in main directory for information on usage and          **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.                  **
**                                                                           **
** Developed by:                                                             **
**      Tan Biao (ginkgoltd@outlook.com)                                     **
**                                                                           **
** ************************************************************************* */

// $Date: 2020-05-11 $
// Written: Tan Biao
// Revised:
//
// Purpose: This file contains the class definition for StructSystem

// The interface:
//

#include "../headers/structsystem.h"
using namespace std;


/*****************************************************************************/
/*                            other functions                                */
/*****************************************************************************/
bool isFileExists(const string &fname)
{
    struct stat buffer;
    if (stat(fname.c_str(), &buffer))
    {
        return false;
    }
    if (!S_ISREG(buffer.st_mode))
    {
        return false;
    }
    return true;
}

bool isDirExists(const string &path)
{
    struct stat buffer;
    if (stat(path.c_str(), &buffer))
    {
        return false;
    }
    if (!S_ISDIR(buffer.st_mode))
    {
        return false;
    }
    return true;
}

bool makedir(const string &path)
{
    mode_t mode = 0755;
    int ret = mkdir(path.c_str(), mode);

    if (ret == 0)
    {
        return true;
    }

    switch (errno)
    {
        case ENOENT:
        {
            // parent didn't exist, try to create it
            {
                int pos = path.find_last_of('/');
                if (pos == std::string::npos)
                {
                    return false;
                }
                if (!makedir(path.substr(0, pos)))
                {
                    return false;
                }
            }
            // now, try to create again
            return 0 == mkdir(path.c_str(), mode);
        }
        case EEXIST:
        {
            return isDirExists(path);
        }
        default:
        {
            return false;
        }
    }
}


/*****************************************************************************/
/*                            helper functions                               */
/*****************************************************************************/
inline bool isFileExist(const string &fname)
{
    struct stat buffer;
    if (stat(fname.c_str(), &buffer))
    {
        return false;
    }
    if (!S_ISREG(buffer.st_mode))
    {
        return false;
    }
    return true;
}

inline bool isDirExist(const string &path)
{
    struct stat buffer;
    if (stat(path.c_str(), &buffer))
    {
        return false;
    }
    if (!S_ISDIR(buffer.st_mode))
    {
        return false;
    }
    return true;
}

inline bool makeDir(const string &path)
{
    mode_t mode = 0755;
    int ret = mkdir(path.c_str(), mode);

    if (ret == 0)
    {
        return true;
    }

    switch (errno)
    {
        case ENOENT:
        {
            // parent didn't exist, try to create it
            {
                int pos = path.find_last_of('/');
                if (pos == std::string::npos)
                {
                    return false;
                }
                if (!makeDir(path.substr(0, pos)))
                {
                    return false;
                }
            }
            // now, try to create again
            return 0 == mkdir(path.c_str(), mode);
        }
        case EEXIST:
        {
            return isDirExist(path);
        }
        default:
        {
            return false;
        }
    }
}

inline void clearCsv(const string &fname)
{
    if (isFileExist(fname))
    {
        fstream fout;
        fout.open(fname, ios::out | ios::trunc);
        fout.close();
    }
}

inline void writeData(fstream &fout, const int ID, const StdArray6d* arr)
{
    fout << ID << ", ";
    for (int i = 0; i < 5; i++)
    {
        fout << (*arr)[i] << ", ";
    }
    fout << (*arr)[5] << endl;
}

double interpolate(vector<double> &xData, vector<double> &yData, double x,
    bool extrapolate)
{
   int size = xData.size();

    // find left end of interval for interpolation
   int i = 0;
   if (x >= xData[size - 2])  // special case: beyond right end
   {
      i = size - 2;
   }
   else
   {
      while (x > xData[i+1]) i++;
   }

    // points on either side (unless beyond ends)
   double xL = xData[i], yL = yData[i], xR = xData[i+1], yR = yData[i+1];
   if (!extrapolate)  // if beyond ends of array and not extrapolating
   {
      if (x < xL) yR = yL;
      if (x > xR) yL = yR;
   }

   double dydx = (yR - yL) / (xR - xL);    // gradient

   return yL + dydx * (x - xL);   // linear interpolation
}


/*****************************************************************************/
/*                               StructSystem                                */
/*****************************************************************************/
StructSystem::StructSystem(int id)
{
    cout << "########## StructSystem Established! ##########" << endl;
    id_ = id;
    h_ = 1e-3;
    alpha_ = 0.2;
    beta_ = 0.0;
    critical_time_step_ = 1e-2;
    workdir_ = getcwd(NULL, 0);
    jobname_ = "filename";
    elem_max_nodes_ = 0;
    elem_types_.push_back("Link");
    elem_types_.push_back("Beam");
    elem_types_.push_back("Plane");
    elem_types_.push_back("Solid");
}

StructSystem::~StructSystem()
{

}

string StructSystem::jobname() const
{
    return jobname_;
}

string StructSystem::workdir() const
{
    return workdir_;
}

double StructSystem::timestep() const
{
    return h_;
}

double StructSystem::critical_time_step() const
{
    return critical_time_step_;
}

const map<int, Particle*>* StructSystem::particles() const
{
    return &particles_;
}

const map<int, BaseMaterial*>* StructSystem::materials() const
{
    return &materials_;
}

const map<int, BaseSection*>* StructSystem::sections() const
{
    return &sections_;
}

const map<int, BaseElement*>* StructSystem::elements() const
{
    return &elements_;
}

const map<int, DOF>* StructSystem::constraints() const
{
    return &constraints_;
}

void StructSystem::info() const
{
    for (int i = 0; i < 80; i++)
    {
        cout << "#";
    }
    cout << endl << "#";
    for (int i = 0; i < 30; i++)
    {
        cout << " ";
    }
    cout << "System Information";
    for (int i = 0; i < 30; i++)
    {
        cout << " ";
    }
    cout << endl;
    for (int i = 0; i < 80; i++)
    {
        cout << "#";
    }
    cout << endl;
    cout << "1. Working Directory: " << workdir_ << endl;
    cout << "2. Job Name: " << jobname_ << endl;
    cout << "3. Number of Particles: " << particles_.size() << endl;
    cout << "4. Number of Elements: " << elements_.size() << endl;
    cout << endl;
    // print paritcle details
    for (int i = 0; i < 31; i++)
    {
        cout << "#";
    }
    cout << " Particle Details ";
    for (int i = 0; i < 31; i++)
    {
        cout << "#";
    }
    cout << endl;
    for (auto iter = particles_.begin(); iter != particles_.end(); ++iter)
    {
        iter->second->info();
    }
    // print element details
    cout << endl;
    for (int i = 0; i < 31; i++)
    {
        cout << "#";
    }
    cout << " Element Details ";
    for (int i = 0; i < 32; i++)
    {
        cout << "#";
    }
    cout << endl;
    for (auto iter = elements_.begin(); iter != elements_.end(); ++iter)
    {
        iter->second->info();
    }
    // print conststaint details
    cout << endl;
    for (int i = 0; i < 30; i++)
    {
        cout << "#";
    }
    cout << " Constraint Details ";
    for (int i = 0; i < 30; i++)
    {
        cout << "#";
    }
    cout << endl;
    cout << "PID |    Ux    |    Uy    |    Uz    |" <<
                 "   Rotx   |   Roty   |   Rotz   |" << endl;
    cout << "    | key, val | key, val | key, val |" <<
                 " key, val | key, val | key, val |" << endl;
    for (auto iter = constraints_.begin(); iter != constraints_.end(); ++iter)
    {
        cout << iter->first << " | ";
        for (int i = 0; i < 6; i++)
        {
            cout << iter->second.key[i] << ", " << iter->second.val[i] << " | ";
        }
        cout << endl;
    }
    cout << endl;
}

void StructSystem::setJobname(const string &job)
{
    jobname_ = job;
}

void StructSystem::setWorkdir(const string &dir)
{
    workdir_ = dir;
}

void StructSystem::setTimeStep(double h)
{
    h_ = h;
}

void StructSystem::setDampCoeff(double alpha, double beta)
{
    alpha_ = alpha;
    beta_ = beta;
}

void StructSystem::addParticle(Particle *p)
{
    int id = p->id();
    if (!particles_.count(id))
    {
        particles_[id] = p;
    }
    else
    {
        cerr << "WARNING: Particle " << id <<
        " already exists and will not add to System" << endl;
    }
}

void StructSystem::addParticles(const vector <Particle *> &ps)
{
    for (int i = 0; i < ps.size(); i++)
    {
        addParticle(ps[i]);
    }
}

void StructSystem::addParticles(const map<int, Particle*> &ps)
{
    for (auto it = ps.begin(); it != ps.end(); it++)
    {
        addParticle(it->second);
    }
}

void StructSystem::addElement(BaseElement* e)
{
    int id = e->id();
    if (!elements_.count(id))
    {
        elements_[id] = e;
        // update critical time step
        if (critical_time_step_ > e->calcCriticalTimeStep())
        {
            critical_time_step_ = e->calcCriticalTimeStep();
        }
        // update max element nodes
        if (elem_max_nodes_ <= e->num_of_particles())
        {
            elem_max_nodes_ = e->num_of_particles();
        }
        // update particle_connected_elems_
        const map<int, Particle*>* ps = e->particles();
        for (auto p = (*ps).begin(); p != (*ps).end(); p++)
        {
           if (particle_connected_elems_.count(p->first))
           {
               particle_connected_elems_[p->first] += 1;
           }
           else
           {
               particle_connected_elems_[p->first] = 1;
           }
        }
    }
    else
    {
        cerr << "WARNING: Element " << id <<
        " already exists and will not add to System" << endl;
    }
}

void StructSystem::addElements(const vector <BaseElement*> &elems)
{
    for (int i = 0; i < elems.size(); i++)
    {
        addElement(elems[i]);
    }
}

void StructSystem::addElements(const map<int, BaseElement*> &elems)
{
    for (auto it = elems.begin(); it != elems.end(); it++)
    {
        addElement(it->second);
    }
}

void StructSystem::addMaterial(BaseMaterial* mat)
{
    int id = mat->id();
    if (!materials_.count(id))
    {
        materials_[id] = mat;
    }
    else
    {
        cerr << "WARNING: Material " << id <<
        " already exists and will not add to System" << endl;
    }
}

void StructSystem::addMaterials(const vector<BaseMaterial*> &mats)
{
    for (int i = 0; i < mats.size(); i++)
    {
        addMaterial(mats[i]);
    }
}

void StructSystem::addMaterials(const map<int, BaseMaterial*> &mats)
{
    for (auto it = mats.begin(); it != mats.end(); it++)
    {
        addMaterial(it->second);
    }
}

void StructSystem::addSection(BaseSection* sect)
{
    int id = sect->id();
    if (!sections_.count(id))
    {
        sections_[id] = sect;
    }
    else
    {
        cerr << "WARNING: Section " << id <<
             " already exists and will not add to System" << endl;
    }
}

void StructSystem::addSections(const vector<BaseSection*> &sects)
{
    for (int i = 0; i < sects.size(); i++)
    {
        addSection(sects[i]);
    }
}

void StructSystem::addSections(const map<int, BaseSection*> &sects)
{
    for (auto it = sects.begin(); it != sects.end(); it++)
    {
        addSection(it->second);
    }
}

void StructSystem::addConstraint(int pid, const DOF &dof)
{
    if (particles_.count(pid))
    {
        if (constraints_.count(pid))
        {
            cerr << "WARNING: " << pid << "already constrainted!" << endl;
        }
        else
        {
            constraints_[pid] = dof;
        }
    }
    else
    {
        cerr << "ERROR: " << pid << "have not found!" << endl;
    }
}

int StructSystem::deleteConstraint(int pid)
{
    int n = 0;
    if (constraints_.count(pid))
    {
        int n = constraints_.erase(pid);
    }
    else
    {
        cout << pid << "have not found in constraints!" << endl;
    }

    return n;
}

void StructSystem::changeConstraint(int pid, const DOF &dof)
{
    if (particles_.count(pid))
    {
        if (constraints_.count(pid))
        {
            constraints_[pid] = dof;
        }
        else
        {
            cout << pid << "have not found in constraints!" << endl;
        }
    }
    else
    {
        cout << pid << "have not found!" << endl;
    }
}

void StructSystem::clearConstraint()
{
    constraints_.erase(constraints_.begin(), constraints_.end());
}

void StructSystem::constraint()
{
    for (auto iter = constraints_.begin(); iter != constraints_.end(); iter++)
    {
        int pid = iter->first;
        particles_[pid]->constraintDof(iter->second);
    }
}

void StructSystem::setParticleDisplace(int pid, const StdArray6d &val)
{
    particles_[pid]->setDisplace(val);
}

void StructSystem::addParticleDisplace(int pid, const StdArray6d &val)
{
    particles_[pid]->addDisplace(val);
}

void StructSystem::setAccelerate(double ax, double ay, double az)
{
    StdArray6d acc = {ax, ay, az, 0, 0, 0};
    StdArray6d force = {0, 0, 0, 0, 0, 0};
    const Mass* m;
    for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
    {
        m = iter->second->mass();
        for (int i = 0; i < 6; i++)
        {
            // Caution: the orient of force and acclerate are opposite
            force[i] = -m->lumped_mass[i] * acc[i];
        }
        iter->second->addForce(force);
    }
}

void StructSystem::setExternalForce(int pid, const StdArray6d &val)
{
    if (particles_.count(pid))
    {
        particles_[pid]->setForce(val);
    }
    else
    {
        cerr << "WARNING: Particle " << pid << " not found!" << endl;
    }
}

void StructSystem::addExternalForce(int pid, const StdArray6d &val)
{
    if (particles_.count(pid))
    {
        particles_[pid]->addForce(val);
    }
    else
    {
        cerr << "WARNING: Particle " << pid << " not found!" << endl;
    }
}

void StructSystem::setInternalForce()
{
    // loop for elements
    for (auto iter = elements_.begin(); iter != elements_.end(); iter++)
    {
        iter->second->setParticleForce();
    }
}

void StructSystem::clearParticleForce()
{
    // clear all force
    for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
    {
        iter->second->clearForce();
    }
}

double StructSystem::autoTimeStep()
{
    h_ = pow(10, floor(log10(critical_time_step_)));
    return h_;
}

void StructSystem::solve(double h, double zeta, bool firstStep)
{
    // loop for particles
    for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
    {
        // cdm(iter->second, h, zeta, firstStep);
        if (firstStep)
        {
            iter->second->solveFirstStep(h, zeta);
        }
        else
        {
            iter->second->solve(h, zeta);
        }
    }
    constraint();
}

void StructSystem::solve(bool firstStep)
{
    // loop for particles
    for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
    {
        // cdm(iter->second, h, zeta, firstStep);
        if (firstStep)
        {
            iter->second->solveFirstStep(h_, alpha_);
        }
        else
        {
            iter->second->solve(h_, alpha_);
        }
    }
    constraint();
}

void StructSystem::saveModel(const string &path)
{
    // if directory does not exist, create it
    if (!isDirExist(path))
    {
        makeDir(path);
    }

    // save particles information
    {
        string fname = path + "/coordinates.csv";
        clearCsv(fname);
        fstream fout;
        fout.open(fname, ios::out | ios::app);
        fout << "PID, x, y, z, thx, thy, thz" << endl;
        for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
        {
            writeData(fout, iter->first, iter->second->coordinate());
        }
        fout.close();
    }

    // save element information
    {
        string fname = path + "/elements.csv";
        clearCsv(fname);
        fstream fout;
        fout.open(fname, ios::out | ios::app);
        fout << "EID, Type, ";
        for (int i = 1; i < elem_max_nodes_; i++)
        {
            fout << "p" << i << ", ";
        }
        fout << "p" << elem_max_nodes_ << endl;
        for (auto iter = elements_.begin(); iter != elements_.end(); iter++)
        {
            const map<int, Particle*>* ps = iter->second->particles();
            fout << iter->first << ", " << iter->second->type() << ", ";
            for (auto p = (*ps).begin(); p != (*ps).end(); p++)
            {
                fout << p->first << ", ";
            }
            fout << "\n";
        }
        fout.close();
    }

    // save constraint information
    {
        string fname = path + "/constraints.csv";
        clearCsv(fname);
        fstream fout;
        fout.open(fname, ios::out | ios::app);
        fout << "PID, Ux, Uy, Uz, Rotx, Roty, Rotz" << endl;

        for (auto dof = constraints_.begin(); dof != constraints_.end(); dof++)
        {
            StdArray6d array {};
            for (int i = 0; i < 6; i++)
            {
                if (dof->second.key[i])
                {
                    array[i] = dof->second.val[i];
                }
                else
                {
                    array[i] = 0.0 / 0.0;
                }
            }
            writeData(fout, dof->first, &array);
        }
        fout.close();
    }
}

void StructSystem::saveParticleResult(const string &path)
{
    // if directory does not exist, create it
    if (!isDirExist(path))
    {
        makeDir(path);
    }

    // create displace.csv
    vector<string> items = {"displace", "velocity", "accelerate"};
    for (int i = 0; i < 3; i++)
    {
        // create file
        string fname = path + "/" + items[i] + ".csv";
        clearCsv(fname);
        fstream fout;
        fout.open(fname, ios::out | ios::app);
        fout << "PID, Ux, Uy, Uz, Rotx, Roty, Rotz" << endl;

        // write data
        for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
        {
            if (i == 0)
            {
                iter->second->outputDisplace(fout);
            }
            else if (i == 1)
            {
                iter->second->outputVelocity(fout);
            }
            else
            {
                iter->second->outputAccelarate(fout);
            }
        }

        // close file
        fout.close();
    }
}

void StructSystem::saveElementResult(const std::string &path)
{
    // if directory does not exist, create it
    if (!isDirExist(path))
    {
        makeDir(path);
    }

    // create link_element_force.csv
    string link_fname = path + "/link_element_force.csv";
    clearCsv(link_fname);
    fstream link_out;
    link_out.open(link_fname, ios::out | ios::app);
    string link_header = "EID, N";
    link_out << link_header << endl;

    // create cable_element_force.csv
    string cable_fname = path + "/cable_element_force.csv";
    clearCsv(cable_fname);
    fstream cable_out;
    cable_out.open(cable_fname, ios::out | ios::app);
    string cable_header = "EID, N";
    cable_out << cable_header << endl;

    // create beam_element_force.csv
    string beam_fname = path + "/beam_element_force.csv";
    clearCsv(beam_fname);
    fstream beam_out;
    beam_out.open(beam_fname, ios::out | ios::app);
    string beam_header = "EID, PID, Fx, SFx, SFy, Mx, My, Mz";
    beam_out << beam_header << endl;

    // write data
    for (auto iter = elements_.begin(); iter != elements_.end(); iter++)
    {
        if (iter->second->type().find("Link") != string::npos)
        {
            iter->second->outputForce(link_out);
        }
        else if (iter->second->type().find("Beam") != string::npos)
        {
            iter->second->outputForce(beam_out);
        }
        else if (iter->second->type().find("Cable") != string::npos)
        {
            iter->second->outputForce(cable_out);
        }
    }

    // close file
    link_out.close();
    beam_out.close();
    cable_out.close();
}

void StructSystem::saveSupportReact(const std::string &path)
{
    // if directory does not exist, create it
    if (!isDirExist(path))
    {
        makeDir(path);
    }

    // create file
    string fname = path + "/" + "support_reaction.csv";
    clearCsv(fname);
    fstream fout;
    fout.open(fname, ios::out | ios::app);

    // write header
    string header = "PID, Fx, Fy, Fz, Mx, My, Mz";
    fout << header << endl;

    // write data
    for (auto iter = constraints_.begin(); iter != constraints_.end(); iter++)
    {
        particles_[iter->first]->outputReactionForce(fout);
    }

    // close file
    fout.close();
}

void StructSystem::saveResult(const string &path)
{
    saveModel(path);
    saveParticleResult(path);
    saveElementResult(path);
    saveSupportReact(path);
}

void StructSystem::releaseContainers()
{
    // release particles_
    for (auto iter = particles_.begin(); iter != particles_.end(); iter++)
    {
        delete iter->second;
        particles_.erase(iter);
    }

    // release elements_
    for (auto iter = elements_.begin(); iter != elements_.end(); iter++)
    {
        delete iter->second;
        elements_.erase(iter);
    }

    // release material
    for (auto iter = materials_.begin(); iter != materials_.end(); iter++)
    {
        delete iter->second;
        materials_.erase(iter);
    }

    // release sections
    for (auto iter = sections_.begin(); iter != sections_.end(); iter++)
    {
        delete iter->second;
        sections_.erase(iter);
    }
}


void StructSystem::killElement(const int id)
{
    if (!elements_.count(id))
    {
        cout << "Warning: Element " << id << " not exists!" << endl;
    }
    else
    {
        // change particle_connected_elems_
        const map<int, Particle*>* ps = elements_[id]->particles();
        for (auto p = (*ps).begin(); p != (*ps).end(); p++)
        {
            if (particle_connected_elems_.count(p->first))
            {
                particle_connected_elems_[p->first] -= 1;
            }
            else
            {
                cout << "Warning: Particle " << p->first << " not found!\n";
            }
        }
        failure_elems_[id] = elements_[id];
        elements_.erase(id);
    }
}

void StructSystem::autoKillElement()
{
    // WARNING: erase data in a loop might cause collapse of program
    for (auto iter = elements_.begin(); iter != elements_.end();)
    {
        if (iter->second->is_failed())
        {
            // change particle_connected_elems_
            const map<int, Particle*>* ps = iter->second->particles();
            for (auto p = (*ps).begin(); p != (*ps).end(); p++)
            {
                if (particle_connected_elems_.count(p->first))
                {
                    particle_connected_elems_[p->first] -= 1;
                }
                else
                {
                    cout << "Warning: Particle " << p->first << " not found!\n";
                }
            }
            failure_elems_[iter->first] = iter->second;
            elements_.erase(iter++);
        }
        else
        {
            iter++;
        }
    }
}

// WARNING: the function is integrated in addElement methond, so it will be
// removed in future
void StructSystem::countParticleConnectedElements()
{
    for (auto iter = elements_.begin(); iter != elements_.end(); iter++)
    {
        const map<int, Particle*>* ps = iter->second->particles();
        for (auto p = (*ps).begin(); p != (*ps).end(); p++)
        {
           if (particle_connected_elems_.count(p->first))
           {
               particle_connected_elems_[p->first] += 1;
           }
           else
           {
               particle_connected_elems_[p->first] = 1;
           }
        }
    }
}

void StructSystem::removeFreeParticle()
{
    for (auto iter = particle_connected_elems_.begin();
        iter != particle_connected_elems_.end(); iter++)
    {
        if ((particle_connected_elems_[iter->first] <= 0) &&
            (!constraints_.count(iter->first)))
        {
            free_particles_[iter->first] = particles_[iter->first];
            particles_.erase(iter->first);
        }
    }
}
